Spin-spin correlation function of the 2D XY model with weak site or bond dilution 



o 

(N 



00 

o 

CD 



o3 



o 

>: 

OV 

On 1 
m 



x 



Oleksandr KapikraniarQ 

Institute for Condensed Matter Physics, National Academy of Sciences of Ukraine, 79011 Lviv, Ukraine 

(Dated: March 9, 2012) 

The spin-spin correlation function of the 2D XY model decays as a power law at all temperatures 
below the Berezinskii-Kosterlitz-Thouless transition point with a temperature dependent exponent 
rj — n(T/J) (J is the ferromagnetic coupling strength). It is known from computer experiments 
that in the 2D XY model with site or bond dilution this exponent depends on concentration p of 
removed sites/bonds as well. Knowing the slope dij/dp at point p — 0, one can predict the value 
of the exponent for small dilution concentrations: rj(p) ~ r)(0) + p(drj/dp)\ p =o. As it is shown in 
this paper, the spin- wave Hamiltonian allows to obtain exact results for this slope: (dn / dp)\ p= o = 
T/(2J) + 0{{T/J) 2 ) and T/(ttJ) + 0{{T/J) 2 ) for site and for bond dilution, respectively. 



I. INTRODUCTION 

An effect produced by introduction of structural ran- 
domness is perhaps one of the first aspects one would be 
willing to investigate, once the properties of the model of 
interest have been successfully studied on regular struc- 
tures. While computer experiment data keep accumu- 
lating for diverse models with structural disorder, this 
problem is often a real challenge to the theory, though. 

We consider the two-dimensional XY model (some- 
times referred to as the planar rotator model), which 
Hamiltonian is traditionally written as 



H 



cos 



-Or') 



(1) 



<r,r'> 



with the sum spanning the pairs of nearest neighbors in 
a square lattice of N sites, J > being the coupling 
strength, and the polar angle 9 r representing the only 
degree of freedom which can be attributed to a spin of 
unit length rotating in a plane. 

The 2D XY model is remarkable for its critical prop- 
erties, as this particular combination of lattice dimen- 
sionality and spin symmetry leads to the existence of 
a finite range of temperatures in which the system 
exhibits critical-like behaviour [Berezinskii-Kosterlitz- 
Thouless (BKT) phase] [l|, most notably, the spin- 
spin correlation function decays as a power law with a 
temperature dependent exponent rj = rj(T/J) below the 
BKT transition point Tbkt- 

In the low temperature limit, where the spin-wave ap- 
proximation (SWA) is applicable, i. e. the cosine in the 
Hamiltonian (1) can be replaced by a quadratic expres- 
sion without affecting the system properties significantly, 
one arrives easily at a power law form of the spin-spin 
correlation function, R - ^ , with an exponent linearly de- 
pendent on temperature 0, 0] : 



?7swa = Tjl-nJ. 



(2) 



It is known, however, that as the temperature increases, 
the real exponent increases non-linearly with tempera- 
ture, so that it assumes the exact value of 1/4 at Tbkt 

Given the two both theoretically and experimentally 
(computer experiment is meant here) acknowledged facts 
that the value of the exponent rj at the BKT transition 
point cannot be changed by structural dilution (see, for 
example, [6]) whereas the value of the BKT transition 
temperature is reduced by dilution and depends on its 
concentration [6-8], one can already make a conclusion 
that the value of rj below the BKT point should depend 
not only on temperature but on dilution concentration as 
well. It is also clear that rj should increase with dilution 
concentration for T < Tbkt- It can be interpreted as 
the increase of effective temperature (decrease of effective 
interaction) due to dilution.. 

A number of works have touched this question, mostly 
using computer simulations. For site dilution, when some 
fraction of sites is excluded from Hamiltonian (JXJ) , see 0- 
[Tll ] , and for bond dilution, when some fraction of bonds 
is removed from (fT]), see @. 

The present study logically continues the theoretical 
part of [To| . making a significant advance [l2j and cov- 
ering both the site and bond dilution cases. The focus 
is on the behavior of the spin-spin correlation function 
and the searched quantity is the dilution concentration p 
dependent exponent rj of the correlation function power- 
law decay. It is natural to assume that the exponent 
rj = r)(T/J,p) is an analytic function with respect to p, 
away from the percolation threshold. Below, p will de- 
note the fraction of removed bonds or sites, depending 
on what dilution type is considered. Thus, rj can be pre- 
sented as a power series 



V(P) - v(fy +p(dv/dp)\p=o 



(3) 
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For small dilution concentrations p, it is enough to know 
the slope (drj / dp)\ p= o to estimate the value of exponent rj 
with good precision. So, in our derivation we drop terms 
that lead to higher order terms in p in ([5]). 

As more simple and transparent from the technical 
point of view case, bond dilution is considered first in 
Section [TTJ where the spin-spin correlation function is 



2 



calculated up to the contributions linear in dilution con- 
centration p and temperature T / J. The analogous but 
more technically involved derivation for the correlation 
function of a system with site dilution can be found in 
Section inil The final results for the exponent of the spin- 
spin correlation function of the systems with site and 
bond dilution are given, respectively, by Eqs. (1391) and 
(HZJ) (see Fig. HD. 



(all pairs (n, ai), . . . , (r i; ai) are different), where (. . .) 
means averaging with respect to disorder configurations, 

H= II E P(iPr,a})..., 
\r,a p r Q =0,l / 

hereafter referred to as configurational averaging. 



II. 2D XY MODEL WITH BOND DILUTION 

In this section the case of bond dilution in the 2D XY 
model is considered. First, in Subsection III AI the bond 
diluted spin-wave Hamiltonian and the procedure of con- 
figurational averaging are defined. Then, in Subsection 
IIIB1 the spin-spin correlation function is calculated up to 
the contributions linear in dilution concentration p and 
temperature. 



A. Bond diluted Hamiltonian and configurational 
averaging 

Hamiltonian ([TJ in the SWA and with bond dilution 
can be written as 

#b.d. = ^E E (^-e r+u J 2 (l-^ Q ), (4) 

r a— x,y 

where u x = (a, 0), u y — (0, a) (a is the lattice spacing), 
and p r . a = 1 if bond (r, u Q ) is removed and otherwise 
(see Fig. [1]). Then, any thermodynamic quantity charac- 
terizing the system will depend on the particular choice 
of configuration {p r ,a} of the discrete variables. 

One is willing to consider here what is often referred to 
as quenched dilution, i. e. when there is a fixed fraction 
p of removed bonds distributed randomly in the system 
and frozen at their position [l3j]. Meaningful physical 
quantities can be obtained averaging them over the con- 
figurations with a fixed fraction of removed bonds p. For 
a large system one might as well allow all configurations, 
ascribing them a probabilistic weight 



P({ P r,a})=l[l(l-p)(l-Pr 



PPr 



(1 -p)S r ,c,(l-P r , a ) p E r , Q Pr„ 



(5) 



meaning that a bond is removed with probability p, which 
will lead to the fact that only realizations with fraction 
J2r Q Pr,a/(2iV) ~ p (2N is the number of bonds in the 
initial lattice) of removed bonds will make essential con- 
tribution to the averaged quantities, when N — > oo. It 
immediately follows that 
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FIG. 1. The occupation number p r , a (a — x,y) takes value 1 
if bond (r, u a ) is removed and otherwise. 

It is convenient to rewrite Hamiltonian (|4]) in the 
Fourier transformed variables 6*k = ^ r e lkr 9 r as 

flb.d. = H + H({p r , a }), H{{p^ a }) ee J2Pr, a H a {r), 



where 



with 



H = -J^7kM-k 



- k = :- } | sin 2 -|- + sin 2 -y- 



is the Hamiltonian of the undiluted system, and 

-i 2 



H a (r) 



_J_ 

~2N 



L k 



— zki 



(l-e 



— ik a 



(7) 
(8) 

(9) 
(10) 



The sums over k in ([8]) and (jlOl) span the 1st Brillouin 
zone. 

The thermodynamic average of some physical quantity 
A can be written as 



(A) = TreAe 



■/Tree' 



-0H h 



(11) 



Pr 



Pi Pr! ,ai 



■Pr 



P 



(6) 



Since 9^ is a complex variable (for k ^ 0): #k = + *^k' 
Tvg above means 

//•OO /"OO 
d9 Hi dd i de t , (12) 

where B/2 stands for a half of the 1st Brillouin zone 
excluding k = (0£ and in the other half are not 
independent due to the relations: #i k = # k and 0f_ k = 
—Ofy. Note, that it is possible to extend the bounds 
of integration in (|12[) to infinity, since the functions that 
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stand after the trace are always rapidly decaying at /3 J — > tained using the Taylor series representations of the ex- 

oo. ponential and (1 + x)~ 1 functions with respect to powers 

The configurationally averaged value of (A) can be ob- of H({p ra }). The equalities in ((6]) easily lead to 

I 



EE 



r,a r ,a' 



E 



+p n 



E--E 



i— 1 ii — 1 i„_2 — 1 

E E ' ' ' E 



ri, ai r n ,a n J i 1 — li 2 —l i n — i — l 



(n)H^ (r 2 ) ■ ■ ■ ffAr-r'- 1 (r n -i)gAr 1 (r w ) 
(i - - « 2 )! • • • (in-2 - i n -iV-in-i ] - 



+ 



(13) 



r 



where [...]' means that the terms having any coinciding 
pairs of indexes, = rj,a; = otj, are excluded from the 
sums enclosed in brackets. This result will be applied in 
the next subsection to calculate the spin-spin correlation 
function. 



expression. One arrives at (see, for example, [4j or [K 
for details) 



yV(R)7?_ k (R)/7 k -> - In -+ const. (18) 

' ' R-^roo 7T a 



B. Spin-spin correlation function of the bond 
diluted 2D XY model 

The spin-spin correlation function of the XY model 
described by Hamiltonian H can be written as 

G(R) = U (J**-*)) = SR Tree Tr J H (14) 



with 



ry k (R) = (e- 4kR - l) /VN. 



(15) 



For the undiluted system, Eq. (|8|). one can write, since 
e -k = e< L and = usin g the notations of (fT2|) . 

G (R) = mr 6 e-^ J ^^^ + ^ 



(16) 



where r; k and ?y k denote the real and imaginary parts of 
?7k(R). It is straightforward to get from the Gaussian 
integration: 



Gq(R) = exp 



i ^j^7 7k (R)r ? _ k (R)/ 7 k 



k#0 



(17) 



here and below sums over k span the entire 1st Brillouin 
zone except the point k = 0. 

To obtain the asymptotic behaviour of (fP7T) at R —> oo 
one should use the fact that rj^fj-k — ^ sin 2 ^ oscillates 
very fast comparing to and, thus, can be replaced 

by its average value 2/N everywhere expect the region 
close to the singularity point k = 0. In this region, re- 
placing in the thermodynamic limit N — > oo the sum 
with an integral and taking the leading terms of the Tay- 
lor expansion of sin 2 ^ and 7 k , one gets an integrable 



It is easy to see that this asymptotic expression leads to 
a power-law decay of the spin-spin correlation function, 
R~ v , with an exponent given by ([2]). 

For a system with bond dilution the spin-spin correla- 
tion function is given by (TTJJ with H = -ffb.d., Eqs. ([7])- 
(1101) . Applying the scheme of configurational averaging 
described in Subsection III Al to the correlation function, 
one is able to collect the resulting series into the following 
expression: 



G(R) = Go(R)|l +^E f (e^k^e-^ 



xGq^R) (e- pHa( - r) ) o -11+ 0(p 2 ) , (19) 



where the terms of higher order in p are dropped and 
(. . .) denotes thermodynamic averaging with Hamilto- 
nian © of the undiluted system: 



(...) = Tr e e-^°.../Tr e e-^ ff °. 



(20) 



Now, using the Taylor series representation of an ex- 
ponential and the results of Appendix |A1 [Eqs. (|A5I) and 
14[) ]. one obtains for H a (r) given by (ITU|) 



i+E 



(2n-l)!! fl 



i (2n)!l V2 



V2 (21) 



(here and below (2m)!! = J]™ i 2 h (2to-1)!! = ]\™ =1 {2i- 
1), to = 1,2,..., and 0!! = 1). 
In a similar way, using (|A3|) . (|A4I) . and the notation 

I a (r) ^ ^= ]T e " Jkr C 1 " e ~ tkaa ) ^k/7k, (22) 
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one arrives at 



i-iy 



(20 J) 
n=i z=o v K ; 



n—l 



(2n-l)!! /ir 2M) 1 
:)!!(2n-2Z)! U/ [' 



(20 



The unity and the term with Z = n in (|2"5|) give \/z [see 
(12 1 [) ] - Changing index n i = n — l and rearranging the 
terms of the infinite series, one has 



aV^H = Go (R)V2e^. (24) 



(2i)\ ^ (2Z)H2« 
The Taylor series representation of (1 — a;) - ™/ 2 , 



(i - xT n/2 = i + E 



i=i 



(2l-2 + n)V. n t 

(21)!! nil ' ( j 



with x = 1/2 and n = 1 and 2i + 1 was used in (|2~T1) and 
(HH) , respectively. 

Now, having (|2"Tj) and (l24l) . one can write the spin-spin 
correlation function in the low temperature limit as 

G(R) = G (R)|l- i | 7 E / "( r )} 



G (R)e-w£, Q ^(r) 



(26) 



Noticing that E r ,J«W = 2 £ k ?7k?/-k/7k, from (HID 
immediately follows a power law decay of the correlation 
function, Rr v , with a dilution concentration dependent 
exponent 

Vb.d.(p) = f?(0) + + 0( P 2 ) + 0((T/J) 2 ) 

~ r,(0) (1 + 2p) , (27) 
where r)(0) is the exponent of the pure system, Eq. @. 



way that the four bonds adjacent to each spinless site 
must be removed, so the occupation number 



Pr = 



( 1, if there is no spin on site r; 
I 0, otherwise, 

has to be introduced; then, 



(28) 



H sA . = H Q + H ({ Pr }), H({ Pr }) = Y,PrHi(r), (29) 

r 

where Hq is the Hamiltonian of the pure model, Eq. ©, 
and 



J 



ffi(r) = - ? 5> r -0 r+u ) 2 



(30) 



with u = (±a,0), (0,±a), which in the Fourier variables 
reads as 



iV 



(31) 



k,k' 



with 5k,k' = 7k,k' - 7k - 7k' [7k was defined in ||SJ)]. 

One can notice that expression (|29[) is not precise when 
there are neighboring spin vacancies; in this case, the 
common bond between the vacant sites is subtracted 
from the "pure" Hamiltonian twice, so it is, in fact, 
brought back with an opposite sign. The precise form 
of H ({p r }) would be 



ff({Pr})=£Prffl(r)+ J2 
r <r,r') 



p rPr ,H 2 (r,r'), (32) 



where H 2 (r,r') = ^(9 r — 8 r <) 2 . However, it is not only 
that the second term in ([32]) gives contributions of order 
of p 2 and higher, after configurational averaging, but it 
can be always dropped when considering the spin-spin 
correlation function, since any non-physical extra bonds 
corresponding to neighboring spinless sites in (|29|) are 
isolated from the rest of the system. 



B. Spin-spin correlation function of the site diluted 
2D XY model 



III. 2D XY MODEL WITH SITE DILUTION 

In this section the case of site dilution in the 2D XY 
model is considered. In Subsection MI A[ the site diluted 
spin-wave Hamiltonian is defined, then, in Subsection 
IIIIB1 the spin-spin correlation function is calculated up 
to the contributions linear in dilution concentration p and 
temperature. 



A. Hamiltonian of the 2D XY model with spin 
vacancies 

The spin- wave Hamiltonian of a system with site dilu- 
tion differs from that of bond dilution, Eq. (Q}, in the 



Now, everything said in Section III Al about the bond 
dilution and configurational averaging can be applied to 
site dilution as well with the only difference that here 
occupation numbers p r are defined for each site r, and 
p = Pr~ ~ ^2 r Pr/N is now the fraction (concentration) of 
removed sites. 

Then, dropping the higher order terms with respect to 
dilution concentration p, the configurationally averaged 
correlation function can be written as 

G(R) = Go (R) 1 1 + p E ^ (e l e~ fiH 1 (r) ^ 

xG^R) (e-^ l(r) )^ - 1 ) + 0(P 2 )](33) 
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with 77k given by (|T5j) . 

The thermodynamic averages in (|33[) can be calcu- 
lated using the Taylor series expansion: e ~P Hl ^ = 
Y^ = q{~~ fiH\{r)) n /n\. Then, the problem reduces to the 
calculation of the quantity (e J ^ k Vk8k H™(r)} with ??k 
given by (fT5|) and 7/k = 0, which is presented in Appendix 
[Bl Looking at the results (IB6|) and (|B8I) , it is easy to see 
that 

0000_.y j \ I CO j 

- mw )„=ni:s <-!>•§ =»pd-')'|. 



i=l Z=0 



and, similarly: 

oo ; 

iEk^e-^W^ = Go (R)exp]T(-l)*| 



x exp 



1 OG 

•— — Vf-ivr 



Explicit expressions for the quantities U and I* are given 
in (|Ul ]) -pi l) . 

Finally, from (|33]>, 



g(r) = g (r; 



(34) 

Using the result of Appendix [Cj Eq. f|C5j) . with 77k given 
by (fl"5j) . one has 

G(R) = G (R)|l-2p + p ^ (Vw^ R ) 

I r^0,R 

xe -W,-2,ftM) ^ 

where 

F,(r,R) = [S l (x~X,y~Y)~S l (x,y)} 2 

+ [S l {y~Y,x-X)-S. l {y 1 x)] 2 (36) 
(i = 1, 2) with the functions Si, 6*2 defined in jDl]), (|D2|) . 



Now, one can expand the exponential function, retain- 
ing only the term linear in 1/(3 J: 

G7R) = Go(R)|l-2p-p £ (^b Fl ( r ' R ) 



8/3 J(tt - 2 



■F 2 (r,R) +0((/?J)- 2 ) . (37) 



Then, using the asymptotic forms (|D4[) and (|D5[) . and 
replacing the sum with an integral, one can show that, 
when R = v 'X 2 + y 2 — > 00, the leading term comes from 
the integral which in polar coordinates reads as 



-i f drFi(r,R) 



I? 2 



rdrdip 



where the integral spans the entire system excluding ar- 
eas close to r = and r = R. This integration can be 
realized as follows: 



pR—a p2tt pa\fN i<2it 

drdip — >• / dr dip + dr dip 

r/0,R Ja JO J R+a JO 

R+a p2ir-a/R. 

dr / dip. 

R-a Ja/R 



There is no difficulty in finding the integrals above, so, 
finally, one arrives at 

G(R) = Go (R) { 1 - 2p - p^j HR/a) } , 

which can be written for small concentrations p and low 
temperatures l/(/3J) as 



with 



G(R) ~ (1 - 2p) 



T 



R 



-TJs.d. 



(38) 



Vs.±(p) = v(o) +Pjj + 0(P 2 ) + 0((T/J) 2 ) 
?7(0)(H- ^p), 



(39) 



^0, R 7" 2 (r 2 + R 2 — 2rR cos ip) 



where 77(0) is the exponent of the pure system given by 
@. The factor (1 — 2p) in (f3"5]) . which appeared naturally 
from the expansion, is the probability to have both sites 
that stand in the pair correlation function occupied with 
spins: (1 — p) 2 — > 1 — 2p. 

p->0 



IV. CONCLUSIONS 

The spin-spin correlation function of the 2D XY model 
decays as a power law at all temperatures below the 
Berezinskii-Kosterlitz-Thouless transition point with a 
temperature dependent exponent 77 = rj(T/J). In the 2D 
XY model with site or bond dilution this exponent de- 
pends on concentration p of removed sites/bonds as well. 
The knowledge of the slope dr//dp at point p = allows 
to predict the value of the exponent for small dilution 
concentrations: rj(p) ~ 77(0) +p(dr]/dp)\ p= o. The analyt- 
ical derivation, performed here in the low-temperature 
limit, led to (drj/dp)\ p= o = 7r?y(0) and 277(0) for site and 
bond dilution, respectively, where 77(0) = T/2irJ is the 
well known result for the model without dilution. These 
results are illustrated in Fig. [5] 

The positive sign of {drj/dp)\ p= Q was well expected, 
since, as it was mentioned in Introduction, dilution can 
be interpreted as the increase of effective temperature. 
One might be tempted to equate the left sides of l|2"7]) and 
(f3"9"| to the universal value of 7?(Tbkt) = 1/4 and identify 
the T in the right side as the corresponding critical tem- 
peratures for site and bond dilution. Unfortunately, such 
an estimate of Tbkt (p) as a function of p would not be 
quantitatively reasonable, since (|2"T|) . (|5t?)) were obtained 
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sitedil. -MC [10] 
site dil. - Eq.(39) 
bond dil. - Eq.(27) 




FIG. 2. Analytical (lines) and Monte Carlo (squares, site dilu- 
tion only) results for the ratios ?? s .d. (p)/??(0) and r?b.d. (p)/f?(0) 
(p is the concentration of missing spins and bonds, respec- 
tively). Concerning the analytical results one is referred to 
Eqs. (|39[) and (|27[) . The Monte Carlo data are borrowed from 
[10t | and come from simulations with Wolff cluster algorithm 
at T/J = 0.08. 



in the spin-wave approximation and do not hold for T 
close to Tbkt(p)- 

It is worth noting that in order to compare the re- 
sults for site and bond dilutions it may be more instruc- 
tive to express the concentration of spinless sites, p = 
(number of empty sites) /(number of all sites), through 
the actual concentration of missing bonds, p' = ((four 
bonds) x (number of empty sites) ) / (number of all bonds) . 
(The latter relation holds, of course, only under the as- 
sumption of low dilution concentration, when the prob- 
ability to have neighboring spinless sites is negligible.) 
Finally, noting that the total number of bonds in the sys- 
tem is two times the number of all sites, we have p — p' /2. 
Then one shall compare the exponent 



%. d .(p')=^(o)(i+(7r/2y) 



(40) 



and (|27p for p' — p, which means that we look at the sys- 
tems with the same number of missing bonds (although 
in the case of site dilution all missing bonds are con- 
nected in unbreakable groups of four). One can notice 
that ?7b.d. > Vs.d. for the same concentration of missing 
bonds, which is well expected, since the disordering effect 
must be stronger for a completely random distribution 
of removed bonds in comparison to the site dilution case 
when removed bonds are connected in groups of four, and 
only these groups are distributed randomly then. 

It also should be mentioned that, in principle, taking 
higher order terms in dilution concentration p in (|13[) . 
one would expect to arrive at the end at the correlation 
function with exponent rj(j>) represented by a series in 
powers of p divergent at the percolation threshold value 
p = Ppcrc. for the square lattice [which is exactly 1/2 



for bond dilution and ~ 0.41 for site dilution (see, for 
example, jlij)]. It is interesting in that it might give an 
exact value for the site percolation threshold which is not 
known yet. However, it might be as well not possible to 
carry out this calculation in an exact way, due to very 
high complexity. 
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Appendix A: Expression for \ 6wx ■ ■ 
Looking at (IT51) . it is easy to see that 



'kl 



where 



d 



8 _ d 



8 



d 



drjt ' d?7_ 



d 
8i)l 



G (R), 
(Al) 



(A2) 



Here and below, (. ..) stands for the thermodynamic 
averaging with the Hamiltonian of the undiluted system, 
see Eq. ®. 

Noting that J^j- = 2<5k,k' (<5k.k' is Kronecker delta) 
and establishing some simple recurrent relations when 
taking sequential derivatives from (1171) . one relatively 
easy arrives at 



(-1) 



n—l 



>< e ri 

comb.(2n,/) u—1 



7k,„ 



2n-2l 
w = l ^ k P 



(2/3J) 2 "-' 
V-k 



, (A3) 



where the sum spans all distinguishable combi- 

comb. (2n,Z) 

nations of I pairs (k^k^), (k i2 ,kj 2 ), ... (k^kjj [com- 
binations which can be obtained from each other by per- 
mutations of the pairs are not distinguished], which can 
be formed using ki, k.2, . . . , k2„. It is instructive to point 
out that 



E 



i = 



(2n) 



comb.(2n,/) 



(2\y(2n~2l)W: 



(A4) 
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Note, that when — 0, (|A3I) gives 



n x 

r 0] 



(2/3 J) 



1 

e n 



ki„ ,-k, 



comb. (2 r, 



u—1 u 



.(A5) 



Appendix B: Calculation of (e i Ek "^H^{t) 



To calculate the quantity 
e l ^ Vk0k H[\r)) = {J/N) n 



E 

k 1: ...,k 2 



-i(ki + ...+k 2 „)r 



X0ki,k a ' ' ' #k 2 „-i,k 2 „ 



(Bl) 



'k 2 „e 



13)) . Each Kronecker delta from (|A3[) deletes one sum- 
mation index k' from the sum in (|Bf I) and "connects" 
two k's belonging either to one g: 

N^ 1 E #k,k'#k,k'/7k = AT" 1 ^2 5k,-k/7k, 

k,k' k 

or to two different g's: 



v ' ^ ]5*,k.9k',*^k,k'/7k = A 1 y^ff*,kg-k,*/7k- 

k 



_1 E-' 

k,k' 



The former will be symbolically represented as \ <// and 
the latter as g — g. In the same vein, jxt) will denote 

A -1 X)kS*,k?7ke~ lkr /7k- Note also ' tnat 5k,k' = ffk'.k- 

Using (|A3[) and the symbolic notations introduced 
above, one can write (IBII) as a sum of terms which are 



products of non-factorizable "blocks" , \g 
(V x 9 x »7)j (»? x .g - .g x 77), etc.: 



n 



e iE k »AR«( r )\ =G (R) 
/ o 

x(\F) Al (\^F) A2 ---( 



y ^ y 



A 2 



comb.(2n,/) 
A, 



x(tj x g x r)) Xl ■ ■ • (77 x g ■ 



ffX^V (B2) 



To each term of the combinatorial sum "^2 



comb. (2n,l) 



, de- 



fined after (|A3[) , corresponds a certain set of integer num- 
bers {A 1) ...,A i ,A^...A*}, Xi,X*i = 0,1,2,.... How- 
ever, there are many terms corresponding to the same 
set {Ai, . . . , A;, Xi , . . . A* }. Determining the number of 
terms (combinations of "connections" ) in (|B2I) which cor- 
respond to any particular set of A's, one can use the A's 
as summation indexes. Using shorter notations 



h = \9 - 9 - ■ ■ • - r > Ij = V X 9 



j?x?7 
(B3) 



[for explicite expressions for Ii, I* the reader is referred 
to (|CJ1|) - (|CJ3[) ] . one arrives at 



(-1) 



n-l 



(2/3J) 



n-l 



I 00 n 00 

nEnE 

i=l \ i= j=l A*=0 



ilE^+EO'- 1 )^-' 

i=i j=i 

\ 1 ,...,\ n ^Ai j\l t*>*1 



x5 [ a; - (n - I a a ;;;;;; a / A: • • ■ ■ • 

(B4) 



[5k, k' was defined after Eq. (l3"Tj) ] one needs the result of where 
the previous appendix for (6k 1 



S(x) = I „' X , ? , and A Al '"' A " is the combi- 

V ' I 0, I f Ai,...,Ai 

natorial factor given by the number of combinations of 
connections in the sum in (|B2[) corresponding to the set 
{Ai, . . . , Aj, A*, . . . A* }. The upper possible values of A's 
are finite, of course, for finite n, but are not important 
(and so can be put equal to oo for simplicity), since the 
first Kronecker delta in (|B4[) assures that altogether one 
has I connections between g's and the second Kronecker 
delta assures that one has (n — I) pairs of ?7's; any re- 
alizations {Ai, . . . , A;, A*, . . . A*} that do not fulfill this 
conditions do not contribute to the sum. 



Factor A 



A*,..., A* 
Ai,...,A, 



can be found from a simple combinato- 



rial analysis: it is given by the number of ways of dividing 
n elements g into Ai and A^ "blocks" of one g, A2 and A£ 
"blocks" of two g's, and so on, which is given by (blocks 
with the same number of g's are not distinguished) 

n!/(A 1 !A 2 !---A / ! 

x(l!) Al (2!) A2 ■••(/!) Al (l!)^(2!)^ ••■(n!) A -) , 

times the number of ways of connecting g's inside every 
"block". Consider a "block" of g kl ,k' } , 9k 2 ,k^, S^.kJ 
and count in how many ways one can interconnect all g's 
in it: \g — g — ... — g' . The answer will be 2 z ~ 1 (i 



which is the number of permutations i\ divided by 2i, 
since a) it is a cyclic structure (so only one ith part of 
all permutations give distinct combinations of intercon- 
nections, others are their repetitions) and b) the com- 
bination of connections is not changed by inversion of 
the g's' order (hence only one half of the permutations 
must be counted), and multiplied by 2*, since every g 
has two k's by which it can connect. The same reasoning 
leads to 2- 7_1 j! possible combinations of connections in- 
side a "block" ijx g — g — . . . — g xr), since it is not cyclic. 



Eventually, 



Ai!(t!)*< 



L\ A*!(j!) a ; 



n 



(B5) 



8 



Therefore, one has 

e^^mir)) =Gn(R)— V l ~ ' 

i oo 1 / t \ n oo / r* \ A* 



f). 



7 1 



(B6) 



When 7?k = 0, 

(n oo \ / n \ 
IIE M E iA *" n 
i=lAj=0/ \i=l y 



and hence 



n oo / 7- \ A* / n 



(B7) 



K i=lA,=0 v 7 \i=l 



Appendix C: Calculation of h and I* 



(B8) 



The sums U and I* , introduced in Appendix [B] (|B3I) , 
can be written as 



^ = ^E f *-i( k '- k )/^ 



(Cl) 



and 



7 * 4E^( k . k ')^ ^ (C2) 
A fr-! 7k 7k' 



k,k' 



(i > 1) with 

Jfkk')-— \^ gk,-kigki,-k 2 ' ' ' gkj-i.-kjgkj.k' 



ki,...,ki 



7k! • • • 7ki 



(C3) 

for £ > 1 and Ia(k, k') = 5k. k'- One can notice the obvi- 
ous recurrent relation 

/ i+ i(k,k') = I^7 i (k ) -k*) 5k%k V 7k . . (C4) 



In the thermodynamic limit, one can replace the sum 
■h over the 1st Brillouin zone with the integral 

,,, a [ n / a dk T r"l a dk„, and then, noticing; that 

a 2 pTv/a f ir/a 

— I dk„ 



TV 



10 



dk 



sin 



^sin 2 %2 +S i n 2 hg. 7T 



and 



2 /-Jr/a 



r/a 



sin 2 ^ cos 2 h± 



= — - I dk 7 



it /a 



sin 2 M sin 2 V 



sin- ^ ! ! 



y sin ^ +sin 2k^ 2 7T 



one can show that 



i^£k,-k'Sk',k»/7k' - (l - ^) gk,- 



(gk ,-k" + gk,k" 



1 1 

2 ~ 7T 



7k 7k" 



^ $>k,k'gk',k"/7k' = ( 1 - ^ ) gk 



- (gk ,-k" + gk,k") + ( 7i ~ - I ":i<"k" 

7T \ Z 7T 



and 



A 

k 

Then, it is easy to see that 

Ji(k,k') = A i 5 k)( _ 1 )i k / + Bi (,g k ,-k' + gk,k') + Ci7k7k' 

with coefficients A, Bi and d obeying the recursive re- 
lations 



A i+1 = 11 A 



Bi+i = --A + ( 1-- )Bi , 

7T \ 7T . 



C l+1 = Q - ~) (A + 273,) - Ci , 
and A = 1, S = 0, C = 0. Thus, 

A i = [l- 2 



EM^ 



1 - - 



i-l—j 



2 

1 

7T 



1 



« = (-i)-(i-i)D-ir(i-i 



Finally, one can obtain expressions for Ii and /* and 
find that 



E(-irz* 

i=l 



5k,k' 



k,k' 



= < 



k,k' 



^i e -i(k+k r i if ^^ k?7 _ k , e -(k + k> =0 . 

(C5) 



x 

00 , otherwise 



Appendix D: Functions Si(A,B) and S 2 (A,B) 

In this appendix one finds the asymptotic form for the 
functions 



sin ^ cos ^ 



sin 2 ±- 

a=x,y ° ul 2 



— smAk x COS Bk y , 

(Dl) 



_,2 fca-a 



sin 2 

a=x,y hm 2 



— COS Ak x COS -Bfc y , 

where the sums span the 1st Brillouin zone. It turns 
out that simple analytic expressions can be obtained, as- 
suming that at least one of the arguments A, B is large. 
Using the integral [15 1 



cos a; 
x 2 + a 2 



-\a\ 



2\a\ 



one can show that 



Si(A -> 00, B) 



. 7T / a . a 

1 dk v e~ A ° sin "^ cos Bk v 



sinh 



(Wn^) 
2sin%^ 



ir / a 



dk y e Aky cos Bky 



and 



Si(A,B -> 00) = 



7T / a 



k r a « 



X cos ■ 



sin Afca; 



dfc y e Bk:c sva.Ak x 



So, 



7T A 2 + B 2 ' 



(D4) 



when at least one of its arguments A, B is sufficiently 
large. 

In a similar way one can show that 



S 2 (A, B) 



B 2 -A 2 



2tt [A 2 + B 2 Y 



(D5) 



(D3) if at least one of its arguments A, B is sufficiently large. 
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